// Author : Zhuowen Tu

#include "mcmcMoveBoundary.h"
#include "mcmcDiffusion.h"
#include "mcmcBlock.h"
//extern CIParsingDoc *pDocView;


/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
McStackPositions::McStackPositions()
{
	vc_stack.SetDimension(12, 1);

	tocheck.SetRowVector(2);
	tocheck(0) = 3;
	tocheck(1) = 4;

	Empty();
}

void McStackPositions::Empty(void)
{
	for (int i=0; i<vc_stack.size(); i++)
	{
		vc_stack(i).SetX(-1);
		vc_stack(i).SetY(-1);
	}

	head = -1;
	tail = -1;
}
	
void McStackPositions::Insert(const McImagePoint &a)
{
	if (head < 0)
	{
//		assert(tail < 0);
		
		head = 0;
		tail = head;

		vc_stack(head) = a;
	}
	else if ((head+1)%vc_stack.size() == tail)
	{
		// full
  		head = (int)((head+1)%vc_stack.size());
		vc_stack(head) = a;
		tail = (int)((tail+1)%vc_stack.size());
	}
	else
	{
		head = (int)((head+1)%vc_stack.size());
		vc_stack(head) = a;
	}
}

bool McStackPositions::IsLoopingAround(const McImagePoint &a)
{
	// push to the stack
	Insert(a);

	for (int i=0; i<tocheck.size(); i++)
	{
		if (vc_stack(head) == vc_stack((head+2*vc_stack.size()-tocheck(i))%vc_stack.size() ) &&
			vc_stack(head) == vc_stack((head+2*vc_stack.size()-2*tocheck(i))%vc_stack.size()))
			return true;
	}

	return false;
}

McImagePoint McStackPositions::PreviousVisit(const int i)
{
	return vc_stack((head-i+2*vc_stack.size()) % vc_stack.size());
}

/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
// construction & destruction
MCMC_Diff_Node::MCMC_Diff_Node()
{
	index = -1;
	id_region = -1;
}

// operators
MCMC_Diff_Node& MCMC_Diff_Node::operator=(const MCMC_Diff_Node &x)
{
	if (this != &x)
	{
		index = x.index;
		id_region = x.id_region;
		id_node = x.id_node;
	}
	return *this;
}

bool MCMC_Diff_Node::operator==(const MCMC_Diff_Node &x)
{
	if (index==x.index && id_region==x.id_region && id_node==x.id_node)
		return true;
	else
		return false;
}

////////////////////////////////////////////
////////////////////////////////////////////
// construction & destruction
MCMC_Diffusion::MCMC_Diffusion()
{
	praster = NULL;
	pmap = NULL;
	Size_Diff_Node = 100;
	pdiff_nodes = new LinkList<MCMC_Diff_Node>[Size_Diff_Node];
	assert(pdiff_nodes != NULL);

	pparticles = NULL;
}

MCMC_Diffusion::~MCMC_Diffusion()
{
	praster = NULL;
	pmap = NULL;
	if (pdiff_nodes != NULL)
		delete []pdiff_nodes;
	pdiff_nodes = NULL;
}

// operatrions
void MCMC_Diffusion::Set(MCMC_Map &map, Raster_Color<uchar> &raster, MCMC_Particle *pparticle1)
{
	pmap = &map;
	praster = &raster;

	// initialization
	assert(pmap != NULL);
	assert(praster != NULL);

	pparticles = pparticle1;

	// corner dynamics

//	corner_dynamics.CornerDetection(raster);
}

void MCMC_Diffusion::SortRegionNodes(LinkList<int> *plist_region)
{
	int				id_region,index,id_node;
	MCMC_AtomicRegion		*pregion;
	MCMC_Diff_Node	diff_node;
	LinkList<int>	list;
	bool			bGo;

	// sort all the nodes on those involved regions
	LayerRegion()->SortRegionNodes(plist_region);

	for (id_node=0; id_node<Size_Diff_Node; id_node++)
		pdiff_nodes[id_node].Empty();

	if (plist_region==NULL)
	{
		// sort all the nodes of each region
		id_region = LayerRegion()->GetFirstAtomicRegion();
		while (id_region >= 0)
		{
			list.Insert(id_region);
			id_region = LayerRegion()->GetNextAtomicRegion();
		}
	}
	else
	{
		bGo = plist_region->GetFirst(id_region);
		while (bGo)
		{
			pregion = LayerRegion()->GetAtomicRegion(id_region);
			if (pregion!=NULL && pregion->Valid())
				list.Insert(id_region);
			bGo = plist_region->GetNext(id_region);
		}
	}

	bGo = list.GetFirst(id_region);
	while (bGo)
	{
		pregion = LayerRegion()->GetAtomicRegion(id_region);
		if(pregion!=NULL && pregion->Valid())
		{
			MCMC_NodeVertex* ptemp;
			index = pregion->GetFirstIndex();
			while (index>=0)
			{
				ptemp = pregion->Get(index);
				assert(ptemp != NULL);

				id_node = ptemp->id_node;
				diff_node.id_region = pregion->Id();
				diff_node.index = index;
				diff_node.id_node = id_node;
				// insert this to the list
				pdiff_nodes[id_node%Size_Diff_Node].Insert(diff_node);

				index = pregion->GetNextIndex();
			}
		}
		bGo = list.GetNext(id_region);
	}
}

int	MCMC_Diffusion::GetSortedIndex(const int id_region,const int id_node)
{
	bool bGo;
	MCMC_Diff_Node diff_node;

	assert(id_node>=0 && id_node<LayerRegion()->MaxNodeSpace());
	bGo = pdiff_nodes[id_node%Size_Diff_Node].GetFirst(diff_node);
	while (bGo)
	{
		if (diff_node.id_region==id_region && diff_node.id_node==id_node)
			return diff_node.index;
		bGo = pdiff_nodes[id_node%Size_Diff_Node].GetNext(diff_node);
	}
	return -1;
}

void MCMC_Diffusion::Competition(const int id_region1)

{
	int						id_node,id_cseg,id_region2,id_ver;
	int						index1,index2;
	MCMC_Vertex				*pver;
	MCMC_Contour_Segment	*pcseg;
	MCMC_AtomicRegion		*pregion1,*pregion2;
	McPoint2D				ptOffset,*pOffsetBuf,*pOffsetBuf1,*pOffsetBuf2;
	McPoint2D				ofStat,ofStat_temp,ofSmooth,ofSmooth_temp;

	// update coordinates
	MCMC_Node				*pnode;
	bool					btextureregion,btextureregion_all;
	RF_uniform				rf_brown;


	// set brownian motion random variable
	if (pmap->GetT()>1000)
		rf_brown.Set(5.0,30.0,1.0);
	else if (pmap->GetT()>100)
		rf_brown.Set(5.0,20.0,1.0);
	else if (pmap->GetT()>10)
		rf_brown.Set(0.0,10.0,1.0);
	else
		rf_brown.Set(0.0,5.0,1.0);

//	rf_brown.Set(0.0,0.0001,0.0001);

	// get the region to be moved first

	pregion1 = LayerRegion()->GetAtomicRegion(id_region1);
	if (pregion1==NULL || !pregion1->Valid())
		return;

	pOffsetBuf = new McPoint2D[pregion1->num_node];
	pOffsetBuf1 = new McPoint2D[pregion1->num_node];
	pOffsetBuf2 = new McPoint2D[pregion1->num_node];

	for (int ii=0; ii<pregion1->num_node; ii++)
	{
		pOffsetBuf[ii].Set(0.0,0.0);
		pOffsetBuf1[ii].Set(0.0,0.0);
		pOffsetBuf2[ii].Set(0.0,0.0);
	}

	// move every node
	for (index1=0; index1<pregion1->num_node; index1++)
	{
		id_node = pregion1->GetSortedNode(index1);
		id_ver = pregion1->GetSortedVertex(index1);

		// initialization
		ofStat.Set(0.0,0.0);
		ofSmooth.Set(0.0,0.0);
		ptOffset.Set(0.0,0.0);

		btextureregion = false;
		btextureregion_all = true;

		if (id_ver<0)
		{
			// get offset set interms of this region
			ptOffset = pregion1->GetDiffusionMove(index1,ofStat,ofSmooth);

			if (pregion1->Model()->Type()==MODEL_TEXTURE || pregion1->Model()->Type()==MODEL_INTENSITY_TEXTURE)
				btextureregion = true;

			if (pregion1->Model()->Type()!=MODEL_TEXTURE && pregion1->Model()->Type()!=MODEL_INTENSITY_TEXTURE)
				btextureregion_all = false;

			int id_node1,id_node2;

			id_node1 = pregion1->GetSortedLastNode(index1);
			id_node2 = pregion1->GetSortedNextNode(index1);


			id_cseg = pregion1->GetSortedSeg(index1);
			pcseg = LayerRegion()->GetCSeg(id_cseg);
			assert(pcseg != NULL);

			id_region2 = pcseg->GetTheOtherRegionId(pregion1->Id());
			if (id_region2>=0)
			{
				pregion2 = LayerRegion()->GetAtomicRegion(id_region2);
				assert(pregion2 != NULL);
				index2 = GetSortedIndex(pregion2->Id(),id_node);

				if (pregion2->Model()->Type()==MODEL_TEXTURE || pregion2->Model()->Type()==MODEL_INTENSITY_TEXTURE)
					btextureregion = true;

				if (pregion2->Model()->Type()!=MODEL_TEXTURE && pregion2->Model()->Type()!=MODEL_INTENSITY_TEXTURE)
					btextureregion_all = false;

				// compute offset

				ptOffset = ptOffset+pregion2->GetDiffusionMove(id_node2,id_node,id_node1,ofStat_temp,ofSmooth_temp);

				ofStat = ofStat+ofStat_temp;
				ofSmooth = ofSmooth+ofSmooth_temp;
			}
			else
			{
			    double dnorm,curv;

				dnorm = pregion1->GetNormal(index1,curv);
				ptOffset.Set(-100.0*cos(-dnorm),-100.0*sin(-dnorm));
			}
		}
		else
		{
			LinkList<int> list_id_region;

			ptOffset.Set(0.0,0.0);

			// this node is a vertex
			pver = LayerRegion()->GetVertex(id_ver);
			assert(pver != NULL);
			id_cseg = pver->GetFirstCSeg();

			while (id_cseg>=0)
			{
				pcseg = LayerRegion()->GetCSeg(id_cseg);
				assert(pcseg != NULL);
				// left region
				id_region2 = pcseg->GetLeftAtomicRegion();

				if (id_region2>=0 && !list_id_region.IsIn(id_region2))
				{
					// insert to the list
					list_id_region.Insert(id_region2);
					// get offset
					pregion2 = LayerRegion()->GetAtomicRegion(id_region2);
					index2 = GetSortedIndex(pregion2->Id(),id_node);
					// this vertex may be the one on the boundary of a hole

					if (index2>=0)
					{
						if (pregion2->Model()->Type()==MODEL_TEXTURE || pregion2->Model()->Type()==MODEL_INTENSITY_TEXTURE)
							btextureregion = true;

						if (pregion2->Model()->Type()!=MODEL_TEXTURE && pregion2->Model()->Type()!=MODEL_INTENSITY_TEXTURE)
							btextureregion_all = false;

						ptOffset = ptOffset+pregion2->GetDiffusionMove(index2,ofStat_temp,ofSmooth_temp);

						ofStat = ofStat+ofStat_temp;
						ofSmooth = ofSmooth+ofSmooth_temp;
					}
				}

				// right region
				id_region2 = pcseg->GetRightAtomicRegion();
				if (id_region2>=0 && !list_id_region.IsIn(id_region2))
				{
					// insert to the list
					list_id_region.Insert(id_region2);
					// get offset

					pregion2 = LayerRegion()->GetAtomicRegion(id_region2);
					index2 = GetSortedIndex(pregion2->Id(),id_node);

					// this vertex may be the one on the boundary of a hole
					if (index2>=0)
					{
						if (pregion2->Model()->Type()==MODEL_TEXTURE || pregion2->Model()->Type()==MODEL_INTENSITY_TEXTURE)
							btextureregion = true;

						ptOffset = ptOffset+pregion2->GetDiffusionMove(index2,ofStat_temp,ofSmooth_temp);
						ofStat = ofStat+ofStat_temp;
						ofSmooth = ofSmooth+ofSmooth_temp;
					}
				}

				// get the next contour segment that connects to this vertex
				id_cseg = pver->GetNextCSeg();
			}
		}

		// in case a texture region is involved
		pOffsetBuf[index1] = ptOffset;
		pOffsetBuf1[index1] = ofStat;
		pOffsetBuf2[index1] = ofSmooth;

		id_node = pregion1->GetSortedNode(index1);
		pnode = LayerRegion()->GetNode(id_node);

		double	dAbs,dmax_off;

		dAbs = sqrt(pow(pOffsetBuf[index1].x,2.0)+pow(pOffsetBuf[index1].y,2.0));
/*		if (pmap->GetT()>1000)
			dmax_off = 1.5;
		else if (pmap->GetT()>100)
			dmax_off = 0.8;
		else if (pmap->GetT()>10)
			dmax_off = 0.4;
		else
			dmax_off = 0.1;
*/

		if (pmap->GetT()>1000)
			dmax_off = 0.3;
		else if (pmap->GetT()>100)
			dmax_off = 0.2;
		else if (pmap->GetT()>10)
			dmax_off = 0.1;
		else
			dmax_off = 0.05;

		if (dAbs>dmax_off)
		{
			pOffsetBuf[index1].x = pOffsetBuf[index1].x/dAbs*dmax_off;
			pOffsetBuf[index1].y = pOffsetBuf[index1].y/dAbs*dmax_off;
		}

		// stay inside
		ptOffset = NewPosition(pnode->GetX(), pnode->GetY(), pOffsetBuf[index1]);
		pregion1->SetOff(index1,ptOffset,pOffsetBuf1[index1],pOffsetBuf2[index1]);
	}

	delete []pOffsetBuf;
	delete []pOffsetBuf1;
	delete []pOffsetBuf2;
}



// stay inside the image
McPoint2D MCMC_Diffusion::NewPosition(const double x_old, const double y_old, McPoint2D &pt_offset)
{
	McPoint2D	pt_offset_real;
	double		x_new,y_new;

	// stay inside
	// x coordinates
	if (x_old<=1.0 || x_old>=(double)pmap->Height()-1.0)
		pt_offset.x = 0.0;

	x_new = x_old+pt_offset.x;
	if (x_new < 1.0)
		x_new = 0.0;
	if (x_new > (double)pmap->Height()-1.0)
		x_new = (double)pmap->Height();

	// y coordinates
	if (y_old<=1.0 || y_old>=(double)pmap->Width()-1.0)
		pt_offset.y = 0.0;

	y_new = y_old+pt_offset.y;
	if (y_new < 1.0)
		y_new = 0.0;

	if (y_new > (double)pmap->Width()-1.0)
		y_new = (double)pmap->Width();

	pt_offset_real.Set(x_new, y_new);

	return pt_offset_real;
}


void MCMC_Diffusion::MoveNodes(const int id_region)
{
	int			id_node,index1,id_ver,id_cseg;
	MCMC_Node	*pnode;
	MCMC_Vertex	*pver;
	McPoint2D	of1,of2,of3;
	MCMC_AtomicRegion	*pregion1;

	cout << "Diffusion MoveNodes" << endl; // Siyuan

	// get the region to be moved first
	pregion1 = LayerRegion()->GetAtomicRegion(id_region);
	if (pregion1==NULL || !pregion1->Valid())
		return;

	// get all the segments in sequence before movement
	list_old_cseg.Empty();
	id_cseg = pregion1->GetFirstCSeg();
	while (id_cseg>=0)
	{
		list_old_cseg.Insert(id_cseg);
		id_cseg = pregion1->GetNextCSeg();
	}

	// save all the contour segments that may possiblly be moved
	list_old_cseg.Empty();
	id_ver = pregion1->GetFirstVer();
	while (id_ver>=0)
	{
		pver = LayerRegion()->GetVertex(id_ver);
		if (pver!=NULL)
		{
			id_cseg = pver->GetFirstCSeg();
			while (id_cseg>=0)
			{
				list_old_cseg.Insert_unique(id_cseg);
				id_cseg = pver->GetNextCSeg();
			}
		}
		id_ver = pregion1->GetNextVer();
	}

	// move every node
	for (index1=0; index1<pregion1->num_node; index1++)
	{
		id_node = pregion1->GetSortedNode(index1);
		pnode = LayerRegion()->GetNode(id_node);

		of1 = pregion1->GetOff(index1,of2,of3);

		pnode->SetX(of1.x);
		pnode->SetY(of1.y);

		// if this node is a vertex
		id_ver = pregion1->GetSortedVertex(index1);
		if (id_ver>=0)
		{
			pver = LayerRegion()->GetVertex(id_ver);
			assert(pver != NULL);
			pver->SetX(pnode->GetX());
			pver->SetY(pnode->GetY());
		}
	}
}

void MCMC_Diffusion::Go(LinkList<int> &list_regions)
{
	int			total,i,forall,k=0;

	cout << "Diffusion Go" << endl; // Siyuan
	//pDocView->OnRepaintAllViews();  // Jiadi Yang

	total = 20;
	forall = 5;


	for (i=0; i<total; i++)
	{
		if (i<forall)
		{
			if (i%5==0)
			{
				// smooth all the boundaries of the region
				LayerRegion()->Smooth();
				// sort all the nodes in each region
				SortRegionNodes();
			}

			// all regions are diffused
			int id_region1;
			id_region1 = LayerRegion()->GetFirstAtomicRegion();
			while (id_region1>=0)
			{
				// compute the new postion
				DiffuseRegionBoundary(id_region1,false);
				// move nodes
				MoveNodes(id_region1);
				id_region1 = LayerRegion()->GetNextAtomicRegion();
			}
		}
		else
		{
			MCMC_AtomicRegion *pregion;
			bool		bgo;
			int			id_region;

			if (k%7==0)
			{
				// smooth all the boundaries of the region
				LayerRegion()->Smooth();
				// sort all the nodes in each region
				SortRegionNodes();
			}

			k++;
			bgo = list_regions.GetFirst(id_region);
			while (bgo)
			{
				pregion = LayerRegion()->GetAtomicRegion(id_region);
				assert(pregion!=NULL);

				// compute the new postion
				DiffuseRegionBoundary(id_region,false);

				MoveNodes(id_region);

				bgo = list_regions.GetNext(id_region);
			}
		}
		// pDocView->OnRepaintAllViews(); // Jiadi Yang
	}

//	corner_dynamics.OneStep(*pgraph);

}


void MCMC_Diffusion::DiffuseRegionBoundary(const int id_region, const bool bResort)
{
	MCMC_AtomicRegion *pregion;

	pregion = LayerRegion()->GetAtomicRegion(id_region);
	assert(pregion != NULL);
	if (pregion->Valid())
	{
		if (bResort)
		{
			int				id_ver,id_cseg,i,id_nb;
			MCMC_Vertex		*pver;
			LinkList<int>	list_region;
			MCMC_Contour_Segment *pcseg;

			// retrive all the regions
			id_ver = pregion->GetFirstVer();
			while (id_ver>=0)
			{
				pver = LayerRegion()->GetVertex(id_ver);
				assert(pver != NULL);
				id_cseg = pver->GetFirstCSeg();
				while (id_cseg>=0)
				{
					pcseg = LayerRegion()->GetCSeg(id_cseg);
					assert(pcseg != NULL);
					for (i=0; i<2; i++)
					{
						if (pcseg->boundry_of[i]>=0)
							list_region.Insert_unique(pcseg->boundry_of[i]);
					}
					id_cseg = pver->GetNextCSeg();
				}
				id_ver = pregion->GetNextVer();
			}

			id_nb = pregion->GetFirstNeighbor();
			while (id_nb>=0)
			{
				list_region.Insert_unique(id_nb);
				id_nb = pregion->GetNextNeighbor();
			}
			// sort all the nodes of each region
			SortRegionNodes(&list_region);
		}
		// region competition
		Competition(id_region);
	}
}

void MCMC_Diffusion::GridGo(LinkList<int> &list_regions)
{
	assert(LayerRegion()->pmb_procedure->CheckConsistencyOfTopologyRecord());
	LayerRegion()->CheckGRegionTopology();

	//cout << "Diffusion GridGo" << endl; // Siyuan

	for (int i=0; i<2; i++)
	{
		bool	bgo;
		int		indexid_region,j;

		j = 0;
		bgo = list_regions.GetFirst(indexid_region);
		while (bgo)
		{
			// compute the new postion
			GridDiffuseRegionBoundary(indexid_region);

			j++;
			bgo = list_regions.GetNext(indexid_region);
		}

		// pDocView->OnRepaintAllViews(); // Jiadi Yang
	}

	// the neighborhood relation
	LayerRegion()->ReGenerateNbRelationships();

	LayerRegion()->CheckGRegionTopology();
}


McPoint2D pt_move_aa;

void MCMC_Diffusion::GridDiffuseRegionBoundary(const int indexid)
{
	MCMC_RegionElement	*pregion;
	int					i,j,id_node1,id_node2,num_nodes;
	double				dforce;
	MCMC_Contour_Segment *pcseg;
	MCMC_Node			 *pnode1,*pnode2;

	McStackPositions	stack_positions;
	McImagePoint		pt_visiting,pt_to_stop;
	bool				blooping,bgo_different_start,bconsistent;
	RF_uniform			rf_step;

	pregion = LayerRegion()->GetRegionFromIndexId(indexid);
	if (pregion == NULL)
		return;
	
	int max_node = LayerRegion()->MaxNodeSpace();

	rf_step.Set(0, 1, 1);
	rf_step.SetCdf();
	
	bgo_different_start = pregion->UseFirstStartCSeg();
	while (bgo_different_start)
	{
		// where to stop
		stack_positions.Empty();
		num_nodes = 0;
		id_node1 = pregion->GetFirstNode();
		while (id_node1 >= 0)
		{
			pnode1 = LayerRegion()->GetNode(id_node1);
			
			pt_visiting.SetRow(pnode1->GetRow());
			pt_visiting.SetCol(pnode1->GetCol());
			
			stack_positions.Insert(pt_visiting);
			num_nodes++;

			if (num_nodes > max_node)
			{
				pregion->FindTheStartCSeg();
				return;
			}

			id_node1 = pregion->GetNextNode();
		}
		pt_to_stop = stack_positions.PreviousVisit(2);

		// start
		stack_positions.Empty();
		id_node1 = pregion->GetFirstNode();
		dforce = pregion->GridGetDiffusionMove(id_node1, id_node2);

		for (i=0; i<num_nodes*3; i++)
		{
			dforce *= 6.0;

			pnode1 = LayerRegion()->GetNode(id_node1);
			pnode2 = LayerRegion()->GetNode(id_node2);

			if(pnode1==NULL || pnode2 == NULL)
			{
				bconsistent = false;
//				assert(false);
				break;
			}
			else
				bconsistent = true;

			pt_move_aa.Set(pnode1->GetX(), pnode1->GetY());

			// current position
			pt_visiting.SetRow(pnode1->GetRow());
			pt_visiting.SetCol(pnode1->GetCol());
			blooping = stack_positions.IsLoopingAround(pt_visiting);
			
			// to the end?
			if (pt_visiting==pt_to_stop && i > 5)
 				break;

			pcseg = LayerRegion()->GetCSeg(LayerRegion()->cseg_on_grid.CSegId(pnode1->GetRow(), pnode1->GetCol(), pnode2->GetRow(), pnode2->GetCol()));
			assert(pcseg != NULL);

			if (dforce >= 1.0)
			{
				if (!LayerRegion()->pmb_procedure->MoveNodes(pcseg, id_node1, id_node2, MB_CSEG_RIGHT))
					id_node1 = id_node2;
				else
					for (j=0; j<4; j++)
					{
						id_node1 = LayerRegion()->pmb_procedure->p_nodes_next_to_move[j];
						if (id_node1 >= 0)
						{
							dforce = pregion->GridGetDiffusionMove(id_node1, id_node2);
							if (id_node2 >= 0)
								break;
						}
					}
			}
			else if (dforce <= -1.0)
			{
				if (!LayerRegion()->pmb_procedure->MoveNodes(pcseg, id_node1, id_node2, MB_CSEG_LEFT))
					id_node1 = id_node2;
				else
					for (j=3; j>=0; j--)
					{
						id_node1 = LayerRegion()->pmb_procedure->p_nodes_next_to_move[j];
						if (id_node1 >= 0)
						{
							dforce = pregion->GridGetDiffusionMove(id_node1, id_node2);
							if (id_node2 >= 0)
								break;
						}
					}
			}
			else
				id_node1 = id_node2;

			if (id_node1 == id_node2)
				dforce = pregion->GridGetDiffusionMove(id_node1, id_node2);
			
			if (blooping)
			{
				// go for one more step
				id_node1 = id_node2;
				dforce = pregion->GridGetDiffusionMove(id_node1, id_node2);
				if ((int)rf_step.Sampling()==1)
				{
					id_node1 = id_node2;
					dforce = pregion->GridGetDiffusionMove(id_node1, id_node2);
				}
			}
		}
		if (!bconsistent)
			break;
		bgo_different_start = pregion->UseNextStartCSeg();
	}
	
	if (!bconsistent)
		LayerRegion()->CreateTopologyRecords();

	// find all the closed contour segments for the region
	pregion->FindTheStartCSeg();
}


